MICB425_Project1
Team 1: Albert Chang (), Alison Fong (33399149), Karen Lau (16524143), Yaqian Luo (59751503), Jessica Ngo (31837131), Peter (Kiet) Truong (36645133)
version April 25, 2018
Project_1
Abstract
To study the diversity and biochemical responses of microbial communities in the context of oxygen minimum zones (OMZs), Saanich Inlet was used as a model ecosystem from which water samples were collected at seven major depths spanning the oxycline. A metagenomic study was conducted in which genomic DNA was extracted from the water samples, PCR amplified, assembled into contiguous sequences, and processed into operational taxonomic units (OTUs) and amplicon sequence variants (ASVs) using mothur and QIIME2 pipelines. Based on OTU and ASV results, we chose to focus in on Cyanobacteria as our taxon of interest. Both mothur and QIIME2 data produced alpha diversity values based on Shannon’s diversity index that suggest a decrease of Cyanobacterial abundance as depth increases. In addition, Cyanobacterial abundance significantly differs across depth and oxygen concentration according to both mothur and QIIME2 data; specifically, abundance significantly decreases at deeper depths and in environments with lower concentrations of oxygen. Within the Cyanobacteria phylum, there are 15 OTUs and 51 ASVs across all samples. Abundance of five OTUs from the mothur pipeline showed significant changes across both depth and oxygen. In contrast, abundance of none of the ASVs from the QIIME2 pipeline showed significant changes across depth, although 17 ASVs showed significant changes across oxygen concentrations. An increase of Cyanobacterial abundance at shallow depths may be explained in part by their ability to absorb red and orange light at the upper boundaries of the water column. This is also supported by the significant change in oxygen and chlorophyll A concentrations across the depth profile. Another explanation for low Cyanobacterial abundance at lower depths may be due to changes in temperature; where temperature drops below 10oC at 100m, Cyanobacterial growth stops completely.To study the diversity and biochemical responses of microbial communities in the context of oxygen minimum zones (OMZs), Saanich Inlet was used as a model ecosystem from which water samples were collected at seven major depths spanning the oxycline. A metagenomic study was conducted in which genomic DNA was extracted from the water samples, PCR amplified, assembled into contiguous sequences, and processed into operational taxonomic units (OTUs) and amplicon sequence variants (ASVs) using mothur and QIIME2 pipelines. Based on OTU and ASV results, we chose to focus in on Cyanobacteria as our taxon of interest. Both mothur and QIIME2 data produced alpha diversity values based on Shannon’s diversity index that suggest a decrease of Cyanobacterial abundance as depth increases. In addition, Cyanobacterial abundance significantly differs across depth and oxygen concentration according to both mothur and QIIME2 data; specifically, abundance significantly decreases at deeper depths and in environments with lower concentrations of oxygen. Within the Cyanobacteria phylum, there are 15 OTUs and 51 ASVs across all samples. Abundance of five OTUs from the mothur pipeline showed significant changes across both depth and oxygen. In contrast, abundance of none of the ASVs from the QIIME2 pipeline showed significant changes across depth, although 17 ASVs showed significant changes across oxygen concentrations. An increase of Cyanobacterial abundance at shallow depths may be explained in part by their ability to absorb red and orange light at the upper boundaries of the water column. This is also supported by the significant change in oxygen and chlorophyll A concentrations across the depth profile. Another explanation for low Cyanobacterial abundance at lower depths may be due to changes in temperature; where temperature drops below 10oC at 100m, Cyanobacterial growth stops completely.
Introduction
Oxygen minimum zones (OMZs) are areas in the ocean where dissolved oxygen concentrations fall below 20 \(\mu\)M (1). Due to temperature increases and other effects caused by global warming, OMZs are expanding at a notable rate. Saanich Inlet, a seasonally anoxic fjord off the coast of British Columbia, is a model ecosystem for studying the diversity and biochemical responses of microbial communities to the hypoxic environments commonly observed in OMZs (1, 2). In particular, Saanich Inlet has been used to model the metabolic coupling and symbiotic interactions that occur in OMZs (3). The inlet undergoes recurring cycles of water column stratification and deep water renewal, rendering it a model ecosystem for studying microbial responses to changes in ocean deoxygenation levels (4). Increased levels of primary productivity in ocean surfaces during the spring season, as well as the limited mixing which occurs between the basin and surface waters both result in the development of an anoxic body of water with increasing depth in the Inlet (2). These anoxic regions become populated with chemolithoautotrophs, and eventually lead to a decrease in aerobically respiring organisms found deeper within these zones. Past studies have demonstrated that these kinds of metabolic shifts generally lead to a loss of nitrogen along with the production of greenhouse gases, most notably methane (CH4) and nitrous oxide (N2O) (1).
In order to investigate the changes that occur in these OMZs, water samples of various depths were collected from Saanich Inlet. Genomic DNA was extracted from these to conduct a metagenomics study, allowing to overcome the barrier of uncultivability of these samples and enable a more thorough exploration of the relationship which exists between the microbes and their communities based on genetic distribution of metabolic processes (5, 6, 7, 8). The extracted DNA is sequenced to generate raw data, which can then be assembled into contiguous sequences. These contigs generated by amplicon sequencing are then compared to a sequence database to determine the microbial taxa present in the environment at each water depth. This involves processing the sequencing data, and there currently exists two methods for this type of data analysis: operational taxonomic units (OTUs) and amplicon sequence variants (ASVs). OTU based pipelines work based on clustering reads which differ by less than a fixed dissimilarity threshold (9). This allows more data to be kept, although some may not be representative of the actual taxa in the community. On the contrary, ASV based pipelines resolve these sequence variants by inferring biological sequences in the sample prior to amplification and possible sequencing errors, and are able to distinguish variants which differ by even one nucleotide (9). This treats each ASV as individual species, though has the potential to discard more data and bias towards sequences that are less error-prone.
The objective of this paper is to analyze the data generated from both OTU and ASV pipelines in order to decide which produces more logical inferences, and ultimately determining which pipeline would be preferred to carry out future analysis of collected water samples. The taxonomy of interest which was selected for this comparison was the phylum Cyanobacteria. Cyanobacteria was selected as there are sufficient numbers of OTUs and ASVs to make sound comparisons, but not so much so that computation-wise it would be infeasible.
Methods
Sampling
Samples were obtained on Saanich Inlet Cruise 72 and taken from seven major depths spanning the oxycline: 10, 100, 120, 135, 150, 165 and 200 m. Waters were filtered, and genomic DNA was extracted. Further sampling details can be found in (3).
DNA Sequencing
Samples were PCR amplified using the 515F and 808R primers, then sequenced according to the standard operating protocol on an Illumina MiSeq platform with Phred33 quality scores.
Data Processing
Sequences were processed using either mothur or QIIME2 as follows:
Mothur Pipeline: mothur was used to clean-up the data. In brief, paired end reads were combined into contigs using their overlapping regions. Low quality sequences, useless sequence data, chimeric sequences and singletons were removed. OTUs were then determined at 97% similarity. OTUs were classified using the SILVA databases, and the taxonomies for each OTU were condensed. The OTU table, taxonomy data and sample metadata were subsequently cleaned up and combined into a phyloseq object.
QIIME2 Pipeline: Demultiplexed sequences were imported into QIIME as manifest reads. QIIME was used to clean up the data along with ASV determination in one step. Sequence quality was visually evaluated, and sequence quality trimming was conducted. All other trimming/filtering parameters were left as default. ASV determination was completed using the Dada2 protocol. ASV classification was completed using the Silva version 119 database at 99% similarity. The ASV table, taxonomy data and sample metadata were subsequently cleaned up and combined into a phyloseq object.
Data Analysis
The aforementioned phyloseq objects were imported into R version 3.4.3 (Windows) or 1.1.383 (Mac). The tidyverse, phyloseq, magrittr, knitr and cowplot packages were loaded and used to complete the data analysis. Data was piped into linear models and ANOVA tests to determine statistical significance at the 95% confidence level.
Environment setup and Data Cleaning
## ── Attaching packages ────────────────────── tidyverse 1.2.1 ──## ✔ ggplot2 2.2.1 ✔ purrr 0.2.4
## ✔ tibble 1.4.2 ✔ dplyr 0.7.4
## ✔ tidyr 0.8.0 ✔ stringr 1.2.0
## ✔ readr 1.1.1 ✔ forcats 0.2.0## ── Conflicts ───────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extract##
## The downloaded binary packages are in
## /var/folders/7c/dzs6m9916vv53_38z2zwxr0h0000gn/T//RtmpmGhgxI/downloaded_packages##
## Attaching package: 'cowplot'## The following object is masked from 'package:ggplot2':
##
## ggsave## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`## ...## 614OTUs were removed because they are no longer
## present in any sample after random subsampling## ...## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`## ...## 6OTUs were removed because they are no longer
## present in any sample after random subsampling## ...Results
- How does microbial community structure change with depth and oxygen concentration?
Alpha-diversity comparison:
Mothur
As shown in Figure 1, alpha-diversity by Shannon’s diversity index based on mothur shows an overall decreasing trend as depth increases. Specifically, a stable trend is observed from depths of 0-100m, starts decreasing at 100-150m and stabilizes at 150-200m. Figure 1 also shows an overall increasing trend as oxygen increases. Specifically, Shannon’s diversity index increases from 2.3-4.25 at oxygen concentrations of 0-40\(\mu\)M and decreases from 4.25-3.9 at oxygen concentrations of 40-220\(\mu\)M. It was also observed that Shannon’s diversity index is higher in oxic conditions (3.84 ? 0.45) than anoxic conditions (2.39 ? 0.07).
QIIME2
The trends of Shannon’s diversity index across depth and oxygen based on QIIME2 data are similar to those of mothur data. However, QIIME2 pipeline produces higher Shannon’s diversity index than mothur does. Shannon’s diversity index increases from 2.9-5.2 at oxygen concentrations of 0-40\(\mu\)M and decreases from 5.2-5.1 at oxygen concentrations of 40-200\(\mu\)M. Figure 2 shows that Shannon’s diversity index is higher in oxic conditions (4.80 ? 0.43) than anoxic conditions (3.15 ? 0.18).
Alpha-diversity of mothur data
## `geom_smooth()` using method = 'loess'
| Statistic | Oxic | Anoxic |
|---|---|---|
| Average | 3.8401008 | 2.3884700 |
| Standard deviation | 0.4523233 | 0.0666717 |
Alpha-diversity of QIIME2 data
## `geom_smooth()` using method = 'loess'
| Statistic | Oxic | Anoxic |
|---|---|---|
| Average | 4.7959464 | 3.1546427 |
| Standard deviation | 0.4296851 | 0.1784873 |
Taxa presence and abundance:
Mothur
31 taxa in the phylum level are detected by mothur pipeline (Figure 3). These taxa, however, have abundance at different magnitudes. Among of them, Proteobacteria is the most predominant phylum in all samples with the highest average abundance over 75. On the contrary, phylum Peregrinlbacteria has abundance no larger than 0.001 in the seven samples (Figure 3). Additionally, different taxa have distinct changes in abundance across depth. For instance, both Thaumarchaeota and Verrucomicrobia reach their maximum abundance at depth of 100m and and gradually decrease in abundance until 200m, while Latescibacteria and Fibrobacteres are almost undetectable in shallow water and their abundances increase dramatically at depth of 200m.
QIIME2
QIIME2 pipeline detects 29 known taxa and unknown taxa in phylum level (Figure 4). Proteobacteria is still the most abundant phylum across samples. Taxa Chlorobi and Candidate division OP3 have the smallest abundance no larger than 0.004. Different pipelines may result in different changes in abundance for the common taxa shared by mothur and QIIME2 data. Although the change patterns of phylum Actinobacteria in abundance across depth are the same in both datasets, Chloroflexi abundance increases gradually with depth in QIIME2 data different from its change pattern in mothur data, in which abundance declines at depth from 100m to 135m and increases gradually until 200m.
- Does your taxon of interest significantly differ in abundance with depth and/or oxygen concentration?
Mothur
The difference in cyanobacteria abundance within depth or oxygen was estimated by the linear model using mothur processed data. The statistical results show that abundance of cyanobacteria is significantly different with depth (p = 0.01263) and oxygen (p = 0.00012). However linear models in Figure 5 indicate completely distinct trends of cyanobacteria abundance across depth and oxygen, where there is a decrease in abundance as depth increases and an increase in abundance as oxygen concentrations increases respectively.
QIIME2
For data processed by QIIME2 pipeline, ANOVA tests indicate cyanobacteria abundance in the seven samples has significantly difference across depth (p = 0.014) and oxygen (p = 0.013).The linear models in Figure 6 reveal that cyanobacteria abundance decreases at deeper water or in the environment with insufficient oxygen.
##
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
##
## Residuals:
## 1 2 6 4 5 3 7
## 46.352 -50.700 20.041 -16.609 -3.284 -39.933 44.132
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 173.5309 39.3958 4.405 0.00699 **
## Depth_m -1.0883 0.2864 -3.800 0.01263 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 42.31 on 5 degrees of freedom
## Multiple R-squared: 0.7428, Adjusted R-squared: 0.6914
## F-statistic: 14.44 on 1 and 5 DF, p-value: 0.01263##
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
##
## Residuals:
## 1 2 6 4 5 3 7
## 6.768 -17.048 19.375 -4.217 12.375 -22.627 5.375
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.3745 7.4690 -0.72 0.504007
## O2_uM 0.9582 0.0885 10.83 0.000117 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.87 on 5 degrees of freedom
## Multiple R-squared: 0.9591, Adjusted R-squared: 0.9509
## F-statistic: 117.2 on 1 and 5 DF, p-value: 0.0001167
##
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
##
## Residuals:
## 1 4 5 3 2 6 7
## 63.270 160.404 72.980 -30.172 -221.273 -44.444 -0.766
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 687.7807 122.7658 5.602 0.0025 **
## Depth_m -3.3051 0.8925 -3.703 0.0140 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 131.8 on 5 degrees of freedom
## Multiple R-squared: 0.7328, Adjusted R-squared: 0.6794
## F-statistic: 13.71 on 1 and 5 DF, p-value: 0.01395##
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
##
## Residuals:
## 1 4 5 3 2 6 7
## 0.5171 190.2269 105.9213 18.5371 -121.0450 -61.0787 -133.0787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 159.0787 57.3201 2.775 0.0391 *
## O2_uM 2.5772 0.6792 3.794 0.0127 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 129.5 on 5 degrees of freedom
## Multiple R-squared: 0.7422, Adjusted R-squared: 0.6907
## F-statistic: 14.4 on 1 and 5 DF, p-value: 0.0127
- Within your taxon, what is the richness (number of OTUs/ASVs)?
Mothur
Across all samples, there are 15 OTUs within cyanobacteria phylum. Table 3 presents the numbers of OTUs within cyanobacteria phylum for each sample. Most of samples contain 3-5 OTUs within cyanobacteria, except Saanich_120 and Saanich_200 which only have 1 and 0 cyanobacteria OTUs respectively.
QIIME2
There are 51 ASVs within cyanobacteria phylum across all samples. The number of ASVs within cyanobacteria phylum for each sample is shown in Table 2. Saanich_010, Saanich 120, and Saanich_135 have relatively high ASV number equal to or over 15. However, it is important to note that there are no singletons within the ASV dataset, and the function used to estimate richness, “estimate_richness” from the phyloseq library is highly dependent on the number of singletons, and warns of unreliable or wrong results in the absence of singletons in the data.
Table showing richness of OTUs using mothur data and ASVs using QIIME2 data
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 15 taxa and 7 samples ]
## sample_data() Sample Data: [ 7 samples by 22 sample variables ]
## tax_table() Taxonomy Table: [ 15 taxa by 7 taxonomic ranks ]## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 51 taxa and 7 samples ]
## sample_data() Sample Data: [ 7 samples by 22 sample variables ]
## tax_table() Taxonomy Table: [ 51 taxa by 7 taxonomic ranks ]## Warning in estimate_richness(., measures = c("Observed")): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.| Depth_m | OTU | ASV | |
|---|---|---|---|
| Saanich_010 | 10 | 5 | 17 |
| Saanich_100 | 100 | 4 | 8 |
| Saanich_120 | 120 | 1 | 15 |
| Saanich_135 | 135 | 4 | 17 |
| Saanich_150 | 150 | 3 | 11 |
| Saanich_165 | 165 | 3 | 5 |
| Saanich_200 | 200 | 0 | 2 |
- Do the abundances of OTUs/ASVs within your taxon of interest change significantly with depth and/or oxygen concentration?
Using the linear model for statistical interpretation, after correcting the p-value for multiple comparisons, the abundance of OTUs 0189, 0658, 1104, 3852, and 4312 from the mothur pipeline within the cyanobacteria phylum changed significantly with both depth and oxygen. Interestingly, the abundance of ASVs in the QIIME2 pipeline did not have any significant changes across the depth profiles after correcting for the p-value. However, there were 17 ASVs that had a significant abundance change with respect to oxygen concentrations.
Mothur
| Estimate | Std. Error | t value | P_value | Adjusted_P | |
|---|---|---|---|---|---|
| Otu0189 | -0.9610475 | 0.2669442 | -3.600181 | 0.0155403 | 0.0491962 |
| Otu0658 | -0.0685106 | 0.0190455 | -3.597207 | 0.0155891 | 0.0491962 |
| Otu1104 | -0.0583306 | 0.0164341 | -3.549356 | 0.0163987 | 0.0491962 |
| Otu3852 | -0.0159083 | 0.0044820 | -3.549356 | 0.0163987 | 0.0491962 |
| Otu4312 | -0.0053028 | 0.0014940 | -3.549356 | 0.0163987 | 0.0491962 |
| Estimate | Std. Error | t value | P_value | Adjusted_P | |
|---|---|---|---|---|---|
| Otu0189 | 0.8586534 | 0.0788686 | 10.88714 | 0.0001136 | 0.0004035 |
| Otu0658 | 0.0611354 | 0.0058159 | 10.51177 | 0.0001345 | 0.0004035 |
| Otu1104 | 0.0522766 | 0.0049096 | 10.64789 | 0.0001264 | 0.0004035 |
| Otu3852 | 0.0142572 | 0.0013390 | 10.64789 | 0.0001264 | 0.0004035 |
| Otu4312 | 0.0047524 | 0.0004463 | 10.64789 | 0.0001264 | 0.0004035 |
QIIME2
## [1] "None of ASV has significantly different abundance acrossing depth with QIIME2 data"
| Estimate | Std. Error | t value | P_value | Adjusted_P | |
|---|---|---|---|---|---|
| Asv12 | 0.0855435 | 0.0080338 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv144 | 0.1568297 | 0.0147287 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv294 | 0.4104610 | 0.0429693 | 9.552421 | 0.0002128 | 0.0006784 |
| Asv404 | 0.1948490 | 0.0182993 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv663 | 0.9749372 | 0.0956990 | 10.187539 | 0.0001564 | 0.0005316 |
| Asv790 | 0.0095048 | 0.0008926 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv945 | 0.0950483 | 0.0089265 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv1055 | 0.0380193 | 0.0035706 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv1085 | 0.1710870 | 0.0160677 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv1141 | 0.0665338 | 0.0062485 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv1209 | 0.0285145 | 0.0026779 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv1454 | 0.1283152 | 0.0120508 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv1578 | 0.0285145 | 0.0026779 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv1728 | 0.2281160 | 0.0214236 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv1817 | 0.2946498 | 0.0276721 | 10.647891 | 0.0001264 | 0.0004605 |
| Asv2018 | 0.3390505 | 0.0720444 | 4.706133 | 0.0053079 | 0.0159238 |
| Asv2336 | 0.1045531 | 0.0098191 | 10.647891 | 0.0001264 | 0.0004605 |
- Are the answers to the above the same using mothur and QIIME2 processed data?
In terms of differences between the two pipelines, the overall trend between both datasets were similar, however the statistical interpretations of the datasets varied. For starters, the Shannon diversity for the whole microbial community processed with mothur was generally smaller than that estimated by QIIME 2 processed data. Additionally, ANOVA tests indicated that the Shannon diversity had no significant change across depth from the mothur pipeline (p = 0.054). In contrast, there was statistical significance with the Shannon diversity across depth from the QIIME2 dataset (p = 0.022). Within the Cyanobacterial taxon, there were 51 ASVs calculated from the QIIME2 pipeline compared to 15 OTUs from the Mothur pipeline. Interestingly, while there is a difference in the richness between both datasets, 5 of the OTUs had a significantly different abundance with respect to depth and oxygen profiles, while none of the ASVs had a significant change in abundance across depth and 17 ASVs which had a significant difference in abundance across oxygen concentrations. Despite all these variations, the Cyanobacterial taxon itself had a significant difference in abundance with respect to depth and oxygen profiles along the water column across both pipelines.
ANOVA on alpha-diversity of mothur data across depth
## Df Sum Sq Mean Sq F value Pr(>F)
## Depth_m 1 2.355 2.3554 6.265 0.0543 .
## Residuals 5 1.880 0.3759
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ANOVA on alpha-diversity of QIIME2 data across depth
## Df Sum Sq Mean Sq F value Pr(>F)
## Depth_m 1 3.563 3.563 10.65 0.0224 *
## Residuals 5 1.673 0.335
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Discussion
The significant differences for Cyanobacteria abundance across depth and oxygen may be caused by photosynthesis and water temperature. Bacteria of the phylum Cyanobacteria obtain their energy through photosynthesis. These phototrophs are characterized by phycocyanin, a bluish pigment, which functions as an auxiliary light-harvesting protein complex to chlorophyll. Phycocyanin absorbs orange and red light at approximately 620nm and fluoresces at about 650nm depending on the species (10). This is particularly interesting because orange and red light are typically absorbed within the first 50m of a water column (11). From our data in Fig 5 & 6, it was found that there are significant differences in the abundance of Cyanobacteria across the Saanich Inlet depth profiles in both the mothur and QIIME2 pipelines. It is likely that this significance exists as a result of Cyanobacteria thriving at the upper boundaries of the water column, where they are still able to absorb red and orange light that is essential for Cyanobacterial photosynthesis. Furthermore, these findings are supported by the oxygen and fluorescence profiles across the water column (Fig 11). There is a significant difference in the concentration of oxygen and chlorophyll A across depth profiles. Higher concentrations of chlorophyll A and oxygen were found within the top 50m, which indicates increased photosynthetic activity. Moreover, there was a significant difference in the abundance of Cyanobacteria across oxygen concentrations and chlorophyll A concentrations for both pipelines; high oxygen and chlorophyll A concentrations were associated with a larger abundance. With Cyanobacterial photosynthesis contributing a substantial proportion of oxygen to Earth’s atmosphere, it is not surprising that a high abundance of Cyanobacteria is associated with a high concentration of oxygen and chlorophyll A, and a shallow depth. Moreover, it has been reported that the growth rate of cyanobacteria is significantly influenced by temperature. Cyanobacteria were observed to have a lower growth rate at colder temperatures. For marine cyanobacteria, the optimal growth temperature range is 20 - 27.5oC; at these temperatures, cyanobacteria grow at a rate of 0.8 d-1 (12). When the temperature dropped to approximately 15oC, the growth rate of cyanobacteria slowed to 0.22 d-1. Interestingly, when temperature dropped below 10oC, cyanobacterial growth rates came to a complete stop (13). Therefore, in our study, the decline of cyanobacteria abundance with depth may at least partly attributed to the decreasing temperature. According temperature data for each sample, the temperature is close to 13oC at 10m, and decreases to about 9oC when at depths below 100m. Hence, the growth rate of cyanobacteria at depths below 100m is substantially slower than the growth rate at the water’s surface. This leads to a lower abundance at lower depths.
Implications of potential differences in pipelines for microbial ecology make it difficult to make conclusive statements in research and discovery, since we become unable to differentiate actual biological differences seen and differences due to a particular pipeline being used. This could also suggest that one pipeline is more suited to the dataset. In fact, this difference could also be exploited, and manipulated so that a pipeline is selected based on the results that it gives, rather than the more appropriate pipeline for the given dataset.
In the context of this project, the main difference between the two pipelines is based on whether the pipeline produces OTUs (mothur) or ASVs (QIIME2). Both of them use different clustering algorithms to determine “true” sequences, and as a result, there are usually far more ASVs in the QIIME2 pipeline than OTUs created with the mothur pipeline. Therefore, when doing downstream analysis of the pipeline results, this may be one of the reasons why there is a large disparity in between the numbers of sequences and quality of the sequences produced even when using the same initial data.
However, when it came to counting abundances within our taxon, cyanobacteria, there was a far lower abundance seen in qiime2 results than mothur results which cannot be explained by having more numerous ASVs than OTUs. Interestingly, when running the estimate_richness function on our qiime2 data for determining abundance, this resulted in the warning that our data provided did not have any singletons (supposedly in the output), and that results are probably unreliable or wrong. This error did not occur with running this function on mothur data. Further analysis of why this error is seen with qiime2 data should be done before fully trusting the results of this function with qiime2 data.
In subsequent analyses with these two pipelines, perhaps a more in-depth analysis of each function available with the phyloseq package should be tested. A dataset of a well-studied and known community could be used so that the output of the functions can be compared with respect to the pipeline used. Comparisons can be made for the outputs of each pipeline, and results from the mothur and QIIME2 pipelines can be assessed with reference to the expected results. Use cases should also be considered, and standard use cases for either pipeline should be indicated so exploitation of pipelines for favourable results does not occur. Mothur may be more suited to a particular dataset, while QIIME2 could be more appropriate for another dataset that should be used with a “denoising” algorithm.
Future directions for this project could possibly involve the analysis of water samples from other OMZs at various depths, along with observation of the Cyanobacteria data present in those samples to see if the relationships are exhibited between oxygen, nitrogen, and the phyla population as with this study. These can also be analyzed once again with both mothur and QIIME2 for further comparison between the pipelines.
References
- Walsh DA, Zaikova E, Howes CG, Song YC, Wright JJ, Tringe SG, Hallam SJ. 2009. Metagenome of a versatile chemolithoautotroph from expanding oceanic dead zones. Science 326(5952): 578-582.
- Torres-Beltr?n M, Hawley AK, Capelle D, Zaikova E, Walsh DA., Mueller A, Finke J. 2017. A compendium of geochemical information from the Saanich Inlet water column. Nature scientific data 4(170159).
- Hawley AK, Torres-Beltr?n M, Zaikova E, Walsh DA, Mueller A, Scofield M, Kheirandish S, Payne C, Pakhomova L, Bhatia M, Shevchek O, Gies EA, Fairley D, Malfatii SA, Norbeck AD, Brewer HM, Pasa-Tolic, L, del Rio TG, Suttle CA, Trige S, Hallam SJ. Data Descriptor: A compendium of multi-omic sequence information from the Saanich Inlet water column. Nature scientific data 4(170160).
- Hallam SJ, Torres-Beltr?n M, Hawley AK. Comment: Monitoring microbial responses to ocean deoxygenation in a model oxygen minimum zone. Nature scientific data 4(170158).
- National Research Council. 2007. The new science of metagenomics: revealing the secrets of our microbial planet. National Academies Press, Washington, DC.
- Bowers RM, Kyrpides NC, Stepanauskas R, Harmon-Smith M, Doud D, Reddy TBK, Tringe SG. 2017. Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nature Biotechnology 35: 725-731.
- Falkowski PG, Fenchel T, Delong EF. 2008. The microbial engines that drive Earth’s biogeochemical cycles. Science 320(5879): 1034-1039.
- Wooley JC, Godzik A, Friedberg I. 2010. A primer on metagenomics. PLoS computational biology 6(2), e1000667.
- Callahan BJ, McMurdie PJ, Holmes SP. 2017. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. The ISME journal 11(12): 2639-2643.
- Simis SG, Huot Y, Babin M, Seppala J, and Metsamaa L. 2012. Optimization of variable fluorescence measurements of phytoplankton communities with cyanobacteria. Photosynthesis research 112(1): 13-30.
- Light in the ocean (n.d.) [Online]. Link: https://manoa.hawaii.edu/exploringourfluidearth/physical/ocean-depths/light-ocean.
- Boyd PW, RynearsonTA, Armstrong EA, Fu F, Hayashi K, Hu Z, Hutchins DA, Kudela RM, Litchman E, Mulholland MR, Passow U, Strzepek RF, Whittaker KA, Yu E, and Thomas MK. (2013). Marine phytoplankton temperature versus growth responses from polar to tropical waters - outcome of a scientific community-wide study. PLoS ONE 8(5): e63091.
- Berg M, Sutula M. 2015. Factors affecting the growth of cyanobacteria with special emphasis on the Sacramento-San Joaquin Delta. Southern California Coastal Water Research Project Technical Report 869.